1 Pre-processing

1.1 Load packages

library(Seurat)
library(SeuratWrappers)
library(Signac)
library(harmony)
library(tidyverse)
library(pheatmap)
library(reshape2)
library(ggpubr)
library(plyr)
library(UCell)
library(qlcMatrix)
library(GOfuncR)
library(Homo.sapiens)
library(biomaRt)
library(viridis)
library(grid)
library(writexl)

set.seed(333)

1.2 Parameters

cell_type = "GCBC"

cols_gcbc <- c(
  "DZ G2M phase" = "#6F66A3",
  "DZ-early S phase" = "#006FC4", 
  "DZ-late S phase" = "#004D89",
  "DZ-non-proliferative" = "#356200", 
  "DZ-LZ transition" = "#D4F0A3",
  "LZ" = "#FFF59A",  
  "LZ-DZ-re-entry commitment" = "#E0EFEF", 
  "LZ-DZ-transition" = "#4EB2FF",
  "LZ-proliferative"  = "#9AD3FF", 
  "MBC" = "#E68098", 
  "prePC" = "#FECFC7")

cels_order <- c("DZ-non-proliferative","DZ-LZ transition","LZ",
                "LZ-DZ-re-entry commitment","LZ-proliferative",
                "LZ-DZ-transition","DZ-early S phase",
                "DZ-late S phase","DZ G2M phase","MBC",
                 "prePC")
# Paths
path_to_obj <- paste0(
  here::here("scATAC-seq/results/R_objects/level_4/"),
  cell_type,
  "/",
  cell_type,
  "_chromVar_CISBP_level_4.rds",
  sep = ""
)

path_to_obj_RNA <- paste0(
  here::here("scRNA-seq/3-clustering/4-level_4/"),
  cell_type,
  "/",
  cell_type,
  "_seu_obj_level_5_eta.rds",
  sep = ""
)

path_to_obj_RNA_DE <- paste0(
  here::here("scRNA-seq/3-clustering/4-level_4/"),
  cell_type,
  "/markers_re_entry.rds",
  sep = "")

path_to_obj_regulons <- paste0(
  here::here("scRNA-seq/3-clustering/4-level_4/"),
  cell_type,
  "/",
  cell_type,
  "_seu_obj_level_5_beta_3P_withpyscenic.rds",
  sep = ""
)

figures_folder <-"/Users/pauli/Desktop/PAPER_FIGURES/GCBC/2.scATAC/"
path_to_regulons <- here::here("scATAC-seq/results/Pyscenic/list_of_regs.rds")
path_to_mat <- here::here("scRNA-seq/results/all_targets_KFKB_mean_expression_per_cluster_eta_noMBC.rds")

1.3 Functions

DARS <- function(ident.1,ident.2){
  DARs <- FindMarkers(
  ident.1 = ident.1,
  ident.2 = ident.2,
  object = seurat,
  min.pct = 0.05,
  test.use = 'LR',
  latent.vars = 'nCount_peaks_redefined')
  
  DARs_filtered <- DARs[DARs$p_val_adj < 0.05,]
return(DARs_filtered)}

DARS_matrix <- function(DARs, group){
  avgexpr_mat <- AverageExpression(
  features = row.names(DARs),
  seurat,
  assays = "peaks_redefined",
  return.seurat = F,
  group.by = group,
  slot = "data")
return(avgexpr_mat)}

DARS_filtered_max <- function(mat, quant, cluster_interest){
    matrix_coefficient <- rowMax(mat) / rowSums(mat)
    regions <- which(matrix_coefficient > quantile(matrix_coefficient, 
                                                   quant))
    filtered_regions <- mat[regions,]
    filtered_regions_interest_clusters <- filtered_regions[which(apply(X = filtered_regions,
        MARGIN = 1,FUN = which.max) %in% cluster_interest),]
    print(nrow(filtered_regions_interest_clusters))
    return(filtered_regions_interest_clusters)}


pheatmap_plot <- function(matrices,n_cutree_rows,
                          annotation_col,mycolors_list){
  colfunc <- colorRampPalette(c("#4575B4", "white","darkgreen"))
  pheatmap(matrices,
      scale = "row",
      annotation_col = annotation_col,
      annotation_colors = mycolors_list,
      color = colfunc(100),
      angle_col = 45,
      show_rownames=F,
      border_color = "white",
      cluster_rows = T,
      cluster_cols = F,
      fontsize_col = 10,
      cutree_rows = n_cutree_rows,
      clustering_distance_cols = "euclidean", 
      clustering_method = "ward.D2")}

2 Load data

2.1 Load scATAC-seq data

seurat <- readRDS(path_to_obj)
seurat
## An object of class Seurat 
## 144021 features across 21447 samples within 2 assays 
## Active assay: peaks_redefined (142257 features, 142254 variable features)
##  1 other assay present: chromvar
##  3 dimensional reductions calculated: umap, lsi, harmony
seurat$level_5 <- revalue(seurat$level_5,
                c("DZ/LZ"="DZ-LZ transition",
                  "DZ_Sphase_HistoneHigh"="DZ-late S phase",
                  "DZ_G2M"="DZ G2M phase",
                  "DZ-nonproliferative"="DZ-non-proliferative",
                  "DZ_Sphase"="DZ-early S phase",
                  "LZ"="LZ",
                  "LZ-DZ-re-entry"="LZ-DZ-re-entry commitment",
                  "LZ-proliferative_BCL2A1neg"="LZ-DZ-transition",
                  "LZ-proliferative_BCL2A1pos"="LZ-proliferative",
                  "MBC-like"="MBC",
                  "PC-precursors"="prePC"))

seurat$level_5 <- factor(seurat$level_5, levels = names(cols_gcbc))

p1 <- DimPlot(seurat, 
    group.by = "level_5",
    cols = cols_gcbc,
    pt.size = 2, raster = T) + NoLegend()

#pdf(paste(figures_folder, "Panel_E.pdf", sep=""),width=4,height=4,paper='special') 
p1 

2.2 Load scRNA-seq data

seurat_RNA <- readRDS(path_to_obj_RNA)
seurat_RNA
## An object of class Seurat 
## 37378 features across 72174 samples within 1 assay 
## Active assay: RNA (37378 features, 0 variable features)
##  3 dimensional reductions calculated: pca, harmony, umap

2.3 Load regulons data

seurat_regulons <- readRDS(path_to_obj_regulons)
seurat_regulons
## An object of class Seurat 
## 37592 features across 61724 samples within 2 assays 
## Active assay: pySCENIC (214 features, 0 variable features)
##  1 other assay present: RNA
##  3 dimensional reductions calculated: pca, harmony, umap

2.3.1 Merging re-entry subtypes

seurat$re_entry_merged <- revalue(seurat$level_5,
          c("DZ-LZ transition"="DZ-LZ transition",
            "DZ-late S phase"="DZ-late S phase",
            "DZ G2M phase"="DZ G2M phase",
            "DZ-non-proliferative"="DZ-non-proliferative",
            "DZ-early S phase"="DZ-early S phase",
            "LZ"="LZ",
            "LZ-DZ-re-entry commitment"="LZ-DZ-re-entry commitment",
            "LZ-DZ-transition"="LZ-DZ-re-entry commitment",
            "LZ-proliferative"="LZ-DZ-re-entry commitment",
            "MBC"="MBC",
            "prePC"="prePC"))

seurat_RNA$merged <- revalue(seurat_RNA$names_level_5_clusters_eta,
            c("DZ_LZ transition"="DZ-LZ transition",
            "DZ cell cycle exit"="DZ-non-proliferative",
            "DZ late Sphase"="DZ-late S phase",
            "G2M early Sphase"="DZ G2M phase",
            "G2M late Sphase"="DZ G2M phase",
            "DZ non proliferative"="DZ-non-proliferative",
            "DZ early Sphase"="DZ-early S phase",
            "LZ"="LZ",
            "LZ_DZ reentry commitment"="LZ-DZ-re-entry commitment",
            "LZ_DZ transition"="LZ-DZ-transition",
            "LZ proliferative"="LZ-proliferative",
            "Precursor MBC"="MBC",
            "Reactivated proliferative MBCs"="MBC",
            "prePC"="prePC"))

seurat_regulons$merged <- revalue(seurat_regulons$names_level_5_clusters_beta,
  c("DZ-cell cycle exit"="DZ-non-proliferative",
    "DZ-nonproliferative"="DZ-non-proliferative",
    "DZ-nonproliferative_FOXP1hi"="DZ-non-proliferative",
    "DZ/LZ"="DZ-LZ transition",
    "LZ"="LZ",
    "LZ-BCL2A1 neg"="LZ",
    "LZ-DZ-re-entry early commitment"="LZ-DZ-re-entry commitment",
    "PC-precursors"="prePC",
    "MBC-like_proli1"="MBC",
    "MBC-like_proli2"="MBC",
    "MBC-like_proli3"="MBC",
    "MBC-like_nonproli"="MBC",
    "MBC-like_FCRL4+"="MBC",
    "LZ-proliferative_BCL2A1pos"="LZ-proliferative",
    "LZ-proliferative_BCL2A1neg"="LZ-DZ-transition",
    "DZ_Sphase_HistoneHigh"="DZ-late S phase",
    "DZ_Sphase"="DZ-early S phase",
    "DZ_G2M_CCNBHigh"="DZ G2M phase",         
    "DZ_G2M_HistoneHigh"="DZ G2M phase"))

3 From the DZ to the LZ

A pairwise differential accessibility analysis between the major biological identities revealed that the DZ-LZ transition was seamless, with not major chromatin changes.

Idents(seurat) <- seurat$re_entry_merged
  
step1 <- DARS(ident.1 = "DZ-non-proliferative", 
              ident.2 = "DZ-LZ transition")

step2 <- DARS(ident.1 = "DZ-LZ transition", 
              ident.2 = "LZ")
step1_mat <- DARS_matrix(DARs = step1,
                         group = "re_entry_merged")
step2_mat <- DARS_matrix(DARs = step2,
                        group = "re_entry_merged")

path_matrix <- rbind(step1_mat$peaks_redefined,
                      step2_mat$peaks_redefined)

3.1 Selecting top candidates regions

path_matrix <- rbind(step1_mat$peaks_redefined,
                      step2_mat$peaks_redefined)

cluster_interest = c(4,5,6)
top_candidates1 <- DARS_filtered_max(path_matrix, 0, 
                                     cluster_interest)
## [1] 3858
top_candidates1_selection <- DARS_filtered_max(top_candidates1[,cluster_interest], 
                            0.6, c(1,2,3))
## [1] 1543

3.2 Visualization of top candidates regions in not merged context

step1_not_merged_mat <- DARS_matrix(step1, 
                                  group ="level_5")

step2_not_merged_mat <- DARS_matrix(step2, 
                                  group ="level_5")

path_not_merged_matrix <- rbind(step1_not_merged_mat$peaks_redefined,step2_not_merged_mat$peaks_redefined)

cell_types <- names(table(seurat$level_5))
annotation_col = data.frame(
                    cell_type = cell_types) 
rownames(annotation_col) <- cell_types

mycolors <- cols_gcbc
names(mycolors) <- cell_types
mycolors_list <- list(mycolors)
names(mycolors_list) <- "cell_type"

#pdf(paste0(figures_folder, "heatmap_DARs_DZtoLZ.pdf"),width=8,height=6,paper='special') 
out_DZ_LZ = pheatmap_plot(path_not_merged_matrix[row.names(top_candidates1_selection),cels_order],
              n_cutree_rows = 3,
              annotation_col, 
              mycolors_list)

3.2.1 Definition of the sub-clusters

cutree_out = cutree(out_DZ_LZ$tree_row, k=3)
cutree_df <- as.data.frame(cbind(names(cutree_out),cutree_out))
cluster1 <- cutree_df[cutree_df$cutree_out == 1,]
cluster2 <- cutree_df[cutree_df$cutree_out == 2,]
cluster3 <- cutree_df[cutree_df$cutree_out == 3,]

nrow(cluster1)
## [1] 668
nrow(cluster2)
## [1] 851
nrow(cluster3)
## [1] 24

4 From the LZ to return to the DZ

A pairwise differential accessibility analysis between the major biological identities revealed an extensive chromatin modulation along the canonical pathway.

Idents(seurat) <- seurat$re_entry_merged

step1_re_entry_merged <- DARS(ident.1 = "LZ", 
                              ident.2 = "LZ-DZ-re-entry commitment")

step2_re_entry_merged  <- DARS(ident.1 = "LZ-DZ-re-entry commitment", 
                               ident.2 = "DZ-early S phase")
step1_re_entry_merged_mat <- DARS_matrix(step1_re_entry_merged, 
                                  group ="re_entry_merged")

step2_re_entry_merged_mat <- DARS_matrix(step2_re_entry_merged, 
                                  group ="re_entry_merged")

path_reentry_merged_matrix <- rbind(step1_re_entry_merged_mat$peaks_redefined,
                             step2_re_entry_merged_mat$peaks_redefined)

4.1 Selecting top candidates regions

cluster_interest = c(6,7,2)

top_candidates_reentry <- DARS_filtered_max(path_reentry_merged_matrix, 0, cluster_interest)
## [1] 2653
top_candidates_reentry_selection <- DARS_filtered_max(top_candidates_reentry[,cluster_interest],0.6, c(1,2,3))
## [1] 1061

4.2 Visualization of top candidates regions in not merged context

step1_re_entry_not_merged_mat <- DARS_matrix(step1_re_entry_merged, 
                                  group ="level_5")

step2_re_entry_not_merged_mat <- DARS_matrix(step2_re_entry_merged, 
                                  group ="level_5")

path_reentry_not_merged_matrix <- rbind(step1_re_entry_not_merged_mat$peaks_redefined,step2_re_entry_not_merged_mat$peaks_redefined)

DZ_LZ_heatmap <- pheatmap_plot(path_reentry_not_merged_matrix[row.names(top_candidates_reentry_selection),cels_order],
              n_cutree_rows = 4,
              annotation_col, 
              mycolors_list)

#pdf(paste0(figures_folder, "heatmap_DARs_LZtoDZ.pdf"),width=8,height=6,paper='special') 
DZ_LZ_heatmap

# Definition of the sub-clusters
cutree_DZ_LZ_heatmap_non_merged = cutree(DZ_LZ_heatmap$tree_row, 
                                         k=4)
cutree_df <- as.data.frame(cbind(names(cutree_DZ_LZ_heatmap_non_merged),cutree_DZ_LZ_heatmap_non_merged))

cluster1_DZ_LZ <- cutree_df[cutree_df$cutree_DZ_LZ_heatmap_non_merged == 1,]
cluster2_DZ_LZ <- cutree_df[cutree_df$cutree_DZ_LZ_heatmap_non_merged == 2,]
cluster3_DZ_LZ <- cutree_df[cutree_df$cutree_DZ_LZ_heatmap_non_merged == 3,]
cluster4_DZ_LZ <- cutree_df[cutree_df$cutree_DZ_LZ_heatmap_non_merged == 4,]

4.2.1 Signature of the LZ cluster

print(paste0("Number of regions:", nrow(cluster2_DZ_LZ)))
## [1] "Number of regions:208"
pheatmap_plot(path_reentry_not_merged_matrix[unlist(cluster2_DZ_LZ$V1),cels_order], n_cutree_rows = 1, annotation_col,  mycolors_list)

seurat <- AddModuleScore(
  seurat,
  name = 'LZ-module',
  features = list(cluster2_DZ_LZ$V1))

colfunc <- colorRampPalette(c("#4575B4", "white","darkgreen"))

#pdf(paste(figures_folder, "Panel_G1.pdf", sep=""),width=5,height=5,paper='special') 
FeaturePlot(seurat,
            order = TRUE,
            min.cutoff = 'q1',
            max.cutoff = 'q99',
            pt.size = 2,
            features  = "LZ.module1",
            raster = T) + 
  theme(plot.title = element_blank()) + 
  scale_color_gradientn( colours = c('lightgrey', 'darkgreen')) + NoLegend()

4.2.1.1 Motif enrichment analysis cluster LZ

enriched.motifs <- FindMotifs(
  object = seurat,
  features = cluster2_DZ_LZ$V1
)

write_xlsx(enriched.motifs,here::here("scATAC-seq/results/tables/supplementary_table_GCBC_MEA_re_entryI.xlsx"))

DT::datatable(enriched.motifs)

4.2.2 Signature cluster re-entry I

print(paste0("Number of regions:", nrow(cluster3_DZ_LZ)))
## [1] "Number of regions:293"
pheatmap_plot(path_reentry_not_merged_matrix[unlist(cluster3_DZ_LZ$V1),cels_order],
              n_cutree_rows = 1,
              annotation_col, 
              mycolors_list)

seurat <- AddModuleScore(
  seurat,
  name = 'reentry-module',
  features = list(cluster3_DZ_LZ$V1))

colfunc <- colorRampPalette(c("#4575B4", "white","darkgreen"))

#pdf(paste(figures_folder, "Panel_G2.pdf", sep=""),width=5,height=5,paper='special') 
FeaturePlot(seurat,
            order = TRUE,
            min.cutoff = 'q1',
            max.cutoff = 'q99',
            pt.size = 2,
            features  = "reentry.module1",
            raster = T) + 
  theme(plot.title = element_blank()) + 
  scale_color_gradientn( colours = c('lightgrey', 'darkgreen')) + NoLegend()

4.2.2.1 Motif enrichment re-entry I

enriched.motifs <- FindMotifs(
  object = seurat,
  features = cluster3_DZ_LZ$V1
)

write_xlsx(enriched.motifs,here::here("scATAC-seq/results/tables/supplementary_table_GCBC_MEA_re_entryII.xlsx"))

DT::datatable(enriched.motifs)

4.2.3 Signature re-entry II

print(paste0("Number of regions:", nrow(cluster1_DZ_LZ)))
## [1] "Number of regions:467"
pheatmap_plot(path_reentry_not_merged_matrix[unlist(cluster1_DZ_LZ$V1),cels_order],
              n_cutree_rows = 1,
              annotation_col, 
              mycolors_list)

seurat <- AddModuleScore(
  seurat,
  name = 'reentryII-module',
  features = list(cluster1_DZ_LZ$V1))

colfunc <- colorRampPalette(c("#4575B4", "white","darkgreen"))

#pdf(paste(figures_folder, "Panel_G3.pdf", sep=""),width=5,height=5,paper='special') 
FeaturePlot(seurat,
            order = TRUE,
            min.cutoff = 'q1',
            max.cutoff = 'q99',
            pt.size = 2,
            features  = "reentryII.module1",
            raster = T) + 
  theme(plot.title = element_blank()) + 
  scale_color_gradientn( colours = c('lightgrey', 'darkgreen')) + NoLegend()

4.2.3.1 Motif enrichment analysis cluster 2

enriched.motifs <- FindMotifs(
  object = seurat,
  features = cluster1_DZ_LZ$V1
)

write_xlsx(enriched.motifs,here::here("scATAC-seq/results/tables/supplementary_table_GCBC_MEA_re_entryIII.xlsx"))

DT::datatable(enriched.motifs)

5 Studying target transcription factors

5.1 BATF family

p1 <- FeaturePlot(seurat,
                  order = TRUE,
                  features = "ENSG00000156127-LINE436-BATF-D-N1",
                  min.cutoff = 'q1',
                  max.cutoff = 'q99',
                  pt.size = 0.4) + 
  theme(plot.title = element_blank()) 

p2 <- FeaturePlot(seurat_RNA,
                  order = TRUE,
                  features = "BATF", 
                  min.cutoff = 'q1',
                  max.cutoff = 'q99',
                  pt.size = 0.4) + 
  theme(plot.title = element_blank()) 

p1 | p2

seurat[["level_5"]] <- factor(unlist(seurat[["level_5"]]), levels= cels_order)
seurat_RNA[["merged"]] <- factor(unlist(seurat_RNA[["merged"]]), levels= cels_order)
seurat_regulons[["merged"]] <- factor(unlist(seurat_regulons[["merged"]]), levels= cels_order)

p1 <- VlnPlot(seurat,
              features = c("ENSG00000156127-LINE435-BATF-D-N1"), 
              group.by = "level_5", pt.size = 0) + NoLegend() + xlab("") + ylab("") +
  theme(legend.position = "none", 
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(), 
        axis.title.y = element_text(size = rel(1), angle = 0), 
        axis.text.y = element_text(size = rel(1))) 

p2 <- VlnPlot(seurat_RNA,
              features = c("BATF"), 
              group.by = "merged", pt.size = 0) + NoLegend() + xlab("") + ylab("") +
  theme(legend.position = "none", 
       # axis.text.x = element_blank(), 
       # axis.ticks.x = element_blank(), 
        axis.title.y = element_text(size = rel(1), angle = 0), 
        axis.text.y = element_text(size = rel(1))) 

#pdf(paste(figures_folder, "Violin_module_final_r-entry.pdf", sep=""),width=4,height=8,paper='special') 
#patchwork::wrap_plots(plotlist = list(p1,p2), ncol = 1)
p1

p2

5.2 NFkB family

p1 <- FeaturePlot(seurat,
                  order = TRUE,
                  features = "ENSG00000109320-LINE3203-NFKB1-D-N1",
                  min.cutoff = 'q1',
                  max.cutoff = 'q99',
                  pt.size = 0.4) + 
  theme(plot.title = element_blank()) 

p2 <- FeaturePlot(seurat_RNA,
                  order = TRUE,
                  features = "NFKB1", 
                  min.cutoff = 'q1',
                  max.cutoff = 'q99',
                  pt.size = 0.4) + 
  theme(plot.title = element_blank()) 

p1 | p2

FeaturePlot(seurat,
            order = TRUE,
            features = c("ENSG00000162924-LINE3217-REL-D-N3",
                         "ENSG00000173039-LINE3227-RELA-D-N13",
                         "ENSG00000104856-LINE3201-RELB-D",
                         "ENSG00000109320-LINE3203-NFKB1-D-N1",
                         "ENSG00000077150-LINE3189-NFKB2-D-N1"),
            min.cutoff = 'q1',
            max.cutoff = 'q99',
            pt.size = 0.4,
            ncol = 2) + 
  theme(plot.title = element_blank()) 

seurat[["level_5"]] <- factor(unlist(seurat[["level_5"]]), levels= cels_order)
seurat_RNA[["merged"]] <- factor(unlist(seurat_RNA[["merged"]]), levels= cels_order)
seurat_regulons[["merged"]] <- factor(unlist(seurat_regulons[["merged"]]), levels= cels_order)

p1 <- VlnPlot(seurat,
              features = c("ENSG00000109320-LINE3203-NFKB1-D-N1"), 
              group.by = "level_5", pt.size = 0) + NoLegend() + xlab("") + ylab("") +
  theme(legend.position = "none", 
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(), 
        axis.title.y = element_text(size = rel(1), angle = 0), 
        axis.text.y = element_text(size = rel(1))) 

p2 <- VlnPlot(seurat_RNA,
              features = c("NFKB1"), 
              group.by = "merged", pt.size = 0) + NoLegend() + xlab("") + ylab("") +
  theme(legend.position = "none", 
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(), 
        axis.title.y = element_text(size = rel(1), angle = 0), 
        axis.text.y = element_text(size = rel(1))) 

p3 <- VlnPlot(seurat_regulons,
              features = c("NFKB1(+)"), 
              group.by = "merged", pt.size = 0) + NoLegend() + xlab("") + ylab("") +
  theme(legend.position = "none", 
        axis.title.y = element_text(size = rel(1), angle = 0), 
        axis.text.y = element_text(size = rel(1))) 

#pdf(paste(figures_folder, "violin_module_r-entry.pdf", sep=""),width=4,height=8,paper='special') 
#patchwork::wrap_plots(plotlist = list(p1,p2,p3), ncol = 1)
p1

p2

p3

6 Session Info

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS/LAPACK: /Users/pauli/opt/anaconda3/envs/Motif_TF/lib/libopenblasp-r0.3.10.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] writexl_1.4.0                           viridis_0.5.1                           viridisLite_0.4.0                       biomaRt_2.44.4                          Homo.sapiens_1.3.1                      TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 org.Hs.eg.db_3.11.4                     GO.db_3.11.4                            OrganismDbi_1.30.0                      GenomicFeatures_1.40.1                  GenomicRanges_1.40.0                    GenomeInfoDb_1.24.2                     AnnotationDbi_1.50.3                    IRanges_2.22.1                          S4Vectors_0.26.0                        Biobase_2.48.0                          BiocGenerics_0.34.0                     GOfuncR_1.8.0                           vioplot_0.3.7                           zoo_1.8-9                               sm_2.2-5.7                              qlcMatrix_0.9.7                         sparsesvd_0.2                           slam_0.1-47                             Matrix_1.3-4                            UCell_1.2.0                             plyr_1.8.6                              ggpubr_0.4.0                            reshape2_1.4.4                         
## [30] pheatmap_1.0.12                         forcats_0.5.0                           stringr_1.4.0                           dplyr_1.0.7                             purrr_0.3.4                             readr_1.4.0                             tidyr_1.1.3                             tibble_3.1.2                            ggplot2_3.3.5                           tidyverse_1.3.0                         harmony_1.0                             Rcpp_1.0.6                              Signac_1.2.1                            SeuratWrappers_0.3.0                    SeuratObject_4.0.2                      Seurat_4.0.3                            BiocStyle_2.16.1                       
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3              SnowballC_0.7.0             rtracklayer_1.48.0          scattermore_0.7             bit64_4.0.5                 knitr_1.30                  irlba_2.3.3                 DelayedArray_0.14.0         data.table_1.14.0           rpart_4.1-15                RCurl_1.98-1.2              generics_0.1.0              cowplot_1.1.1               RSQLite_2.2.1               RANN_2.6.1                  future_1.21.0               bit_4.0.4                   spatstat.data_2.1-0         xml2_1.3.2                  lubridate_1.7.9             httpuv_1.6.1                SummarizedExperiment_1.18.1 assertthat_0.2.1            xfun_0.18                   hms_0.5.3                   evaluate_0.14               promises_1.2.0.1            fansi_0.5.0                 progress_1.2.2              dbplyr_1.4.4                readxl_1.3.1                igraph_1.2.6                DBI_1.1.0                   htmlwidgets_1.5.3           spatstat.geom_2.2-0         ellipsis_0.3.2              crosstalk_1.1.1             backports_1.1.10            bookdown_0.21               deldir_0.2-10               vctrs_0.3.8                 here_1.0.1                 
##  [43] remotes_2.4.1               ROCR_1.0-11                 abind_1.4-5                 withr_2.4.2                 ggforce_0.3.2               sctransform_0.3.2           GenomicAlignments_1.24.0    prettyunits_1.1.1           goftest_1.2-2               cluster_2.1.0               lazyeval_0.2.2              crayon_1.4.1                labeling_0.4.2              pkgconfig_2.0.3             tweenr_1.0.1                nlme_3.1-150                rlang_0.4.11                globals_0.14.0              lifecycle_1.0.0             miniUI_0.1.1.1              BiocFileCache_1.12.1        modelr_0.1.8                rsvd_1.0.3                  rprojroot_2.0.2             cellranger_1.1.0            polyclip_1.10-0             matrixStats_0.59.0          lmtest_0.9-38               graph_1.66.0                ggseqlogo_0.1               carData_3.0-4               reprex_0.3.0                ggridges_0.5.3              png_0.1-7                   bitops_1.0-7                KernSmooth_2.23-17          Biostrings_2.56.0           blob_1.2.1                  parallelly_1.26.1           mapplots_1.5.1              rstatix_0.6.0               ggsignif_0.6.0             
##  [85] scales_1.1.1                memoise_1.1.0               magrittr_2.0.1              ica_1.0-2                   zlibbioc_1.34.0             compiler_4.0.3              RColorBrewer_1.1-2          fitdistrplus_1.1-5          Rsamtools_2.4.0             cli_3.0.0                   XVector_0.28.0              listenv_0.8.0               patchwork_1.1.1             pbapply_1.4-3               MASS_7.3-53                 mgcv_1.8-33                 tidyselect_1.1.1            stringi_1.6.2               yaml_2.2.1                  askpass_1.1                 ggrepel_0.9.1               fastmatch_1.1-0             tools_4.0.3                 future.apply_1.7.0          rio_0.5.16                  rstudioapi_0.11             foreign_0.8-80              lsa_0.73.2                  gridExtra_2.3               farver_2.1.0                Rtsne_0.15                  digest_0.6.27               BiocManager_1.30.10         shiny_1.6.0                 car_3.0-10                  broom_0.7.2                 later_1.2.0                 RcppAnnoy_0.0.18            httr_1.4.2                  colorspace_2.0-2            rvest_0.3.6                 XML_3.99-0.3               
## [127] fs_1.5.0                    tensor_1.5                  reticulate_1.20             splines_4.0.3               uwot_0.1.10                 RBGL_1.64.0                 RcppRoll_0.3.0              spatstat.utils_2.2-0        plotly_4.9.4.1              xtable_1.8-4                jsonlite_1.7.2              R6_2.5.0                    pillar_1.6.1                htmltools_0.5.1.1           mime_0.11                   DT_0.16                     glue_1.4.2                  fastmap_1.1.0               BiocParallel_1.22.0         codetools_0.2-17            utf8_1.2.1                  lattice_0.20-41             spatstat.sparse_2.0-0       curl_4.3.2                  leiden_0.3.8                gtools_3.9.2                zip_2.1.1                   openxlsx_4.2.4              openssl_1.4.4               survival_3.2-7              rmarkdown_2.5               docopt_0.7.1                munsell_0.5.0               GenomeInfoDbData_1.2.3      haven_2.3.1                 gtable_0.3.0                spatstat.core_2.2-0